#pragma once

void compute_delta( O55 D, O55 Delta)
{ 
   int iflag;


        //1. compute Delta 

        // a) LL
          int Nmix=4;
          int Nad=1;
          int LLmixprime[4]={2,3,4,5}; // En c empezamos en 0 asi que habra que restar
          int LLadprime[1]={1};
          double **Dtilde=malloc(sizeof(double*)*(Nmix));
          for (int iop1=0; iop1<Nmix; iop1++) Dtilde[iop1]=malloc(sizeof(double)*(Nmix));
          double **Dtilde_inv=malloc(sizeof(double*)*(Nmix));
          for (int iop1=0; iop1<Nmix; iop1++) Dtilde_inv[iop1]=malloc(sizeof(double)*(Nmix));
          double *Delta_LL=malloc(sizeof(double)*(Nmix));
          for (int i=0; i<Nmix; i++) Delta_LL[i]=0;

        for (int imix=0; imix<Nmix; imix++)
        {
                int imixprime=LLmixprime[imix]-1;
                for (int kmix=0; kmix<Nmix; kmix++)
                {
                        int kmixprime=LLmixprime[kmix]-1;
                        Dtilde[imix][kmix]=D[kmixprime][imixprime];
                }

        }
        inverse_real(Dtilde,Nmix,Dtilde_inv,iflag);
	if (iflag==1) {printf("D matrix is singular\n");}
        for (int iad=0; iad<Nad; iad++)
        {
                int iadprime=LLadprime[iad]-1;

                for(int kmix=0; kmix<Nmix; kmix++)
                {
                        for (int jmix=0; jmix<Nmix; jmix++)
                        {
                                int jmixprime=LLmixprime[jmix]-1;
                                Delta_LL[kmix]=Delta_LL[kmix]-Dtilde_inv[kmix][jmix]*D[iadprime][jmixprime];
                        }
                }
        }
        // b) LR
        Nad=2;
        int LRmixprime[3]={1,4,5};
        int LRadprime[2]={2,3};
        double **Dtildem1=malloc(sizeof(double*)*(Nmix-1));
        for (int iop1=0; iop1<Nmix; iop1++) Dtildem1[iop1]=malloc(sizeof(double)*(Nmix-1));
        double **Dtildem1_inv=malloc(sizeof(double*)*(Nmix-1));
        for (int iop1=0; iop1<Nmix; iop1++) Dtildem1_inv[iop1]=malloc(sizeof(double)*(Nmix-1));
        double **Delta_LR=malloc(sizeof(double*)*(Nad));
        for (int iad=0; iad<Nad; iad++) Delta_LR[iad]=malloc(sizeof(double)*(Nad));
         for (int iad1=0; iad1<Nad; iad1++)  for (int iad2=0; iad2<Nad; iad2++) Delta_LR[iad1][iad2]=0;

        for (int imix=0; imix<Nmix-1; imix++)
        {
                int imixprime=LRmixprime[imix]-1;
                for (int kmix=0; kmix<Nmix-1; kmix++)
                {
                        int kmixprime=LRmixprime[kmix]-1;
                        Dtildem1[imix][kmix]=D[kmixprime][imixprime];
                }
        }
        inverse_real(Dtildem1,Nmix-1,Dtildem1_inv,iflag);
	if (iflag==1) {printf("D matrix is singular\n");}
        for (int iad=0; iad<Nad; iad++)
        {
                int iadprime=LRadprime[iad]-1;
                for (int kmix=0; kmix<Nmix-1; kmix++)
                {
                        for (int jmix=0; jmix<Nmix-1; jmix++)
                        {
                                int jmixprime=LRmixprime[jmix]-1;
                                Delta_LR[iad][kmix]=Delta_LR[iad][kmix]-Dtildem1_inv[kmix][jmix]*D[iadprime][jmixprime];
                        }
                }
        }

        // c) TT
        Nad=2;
        int TTmixprime[3]={1,2,3};
        int TTadprime[2]={4,5};
        double **Delta_TT=malloc(sizeof(double*)*(Nad));
        for (int iad=0; iad<Nad; iad++) Delta_TT[iad]=malloc(sizeof(double)*(Nad));
        for (int iad1=0; iad1<Nad; iad1++)  for (int iad2=0; iad2<Nad; iad2++) Delta_TT[iad1][iad2]=0;

        for (int imix=0; imix<Nmix-1; imix++)
        {
                int imixprime=TTmixprime[imix]-1;
                for (int kmix=0; kmix<Nmix-1; kmix++)
                {
                        int kmixprime=TTmixprime[kmix]-1;
                        Dtildem1[imix][kmix]=D[kmixprime][imixprime];
                }
        }
        inverse_real(Dtildem1,Nmix-1,Dtildem1_inv,iflag);
	if (iflag==1) {printf("D matrix is singular\n");}
        for (int iad=0; iad<Nad; iad++)
        {
                int iadprime=TTadprime[iad]-1;
                for (int kmix=0; kmix<Nmix-1; kmix++)
                {
                        for (int jmix=0; jmix<Nmix-1; jmix++)
                        {
                                int jmixprime=TTmixprime[jmix]-1;
                                Delta_TT[iad][kmix]=Delta_TT[iad][kmix]-Dtildem1_inv[kmix][jmix]*D[iadprime][jmixprime];
                        }
                }
        }



        // 2. Write Delta in the usual matrix form
        Delta[0][0]=0;
        for (int kmix=0; kmix<Nmix; kmix++) { int kmixprime=LLmixprime[kmix]-1; Delta[0][kmixprime]=Delta_LL[kmix];}
        Delta[1][1]=Delta[1][2]=Delta[2][1]=Delta[2][2]=0;
        for (int iad=0; iad<Nad; iad++) for (int kmix=0; kmix<Nmix-1; kmix++)
        {
                int iadprime=LRadprime[iad]-1;
                int kmixprime=LRmixprime[kmix]-1;
                Delta[iadprime][kmixprime]=Delta_LR[iad][kmix];
        }
        Delta[3][3]=Delta[3][4]=Delta[4][3]=Delta[4][4]=0;
        for (int iad=0; iad<Nad; iad++) for (int kmix=0; kmix<Nmix-1; kmix++)
        {
                int iadprime=TTadprime[iad]-1;
                int kmixprime=TTmixprime[kmix]-1;
                Delta[iadprime][kmixprime]=Delta_TT[iad][kmix];
        }

/*        printf("Delta_LL\n");
        for (int kmix=0; kmix<Nmix; kmix++) printf(" %8.5f ", Delta_LL[kmix]);
        printf("\n");
        printf("Delta_LR\n");
         for (int iad=0; iad<Nad; iad++) for (int kmix=0; kmix<Nmix-1; kmix++) printf(" %8.5f ", Delta_LR[iad][kmix]);
        printf("\n");
        printf("Delta_TT\n");
         for (int iad=0; iad<Nad; iad++) for (int kmix=0; kmix<Nmix-1; kmix++) printf(" %8.5f ", Delta_TT[iad][kmix]);
        printf("\n");

        printf("Delta=\n");
        for (int iop1=0; iop1<5; iop1++) {
                for (int iop2=0; iop2<5; iop2++)                printf(" %8.5f ", creal(Delta[iop1][iop2]));
                printf("\n");
*/
}

/*
//Programa eliminacion Gauss   con pivotado: matriz cuadrada y descomposicion LU
void gaussp(int n, double** m, double* v, double** l){
   
//Definir vector  permutacion    
int* per;
per=new int[n];
for (int i=1;i<n;i++) per[i]=i;  


for (int i=0;i<n; i++){
  l[i][i]=1;
  if (i<n-1) for(int j=i+1; j<n; j++) l[i][j]=0;
  }


//ciclo eliminacion  
for (int i=0;i<n;i++){
//search for the pivot
   double piv=abs(m[i][i]);;
   int ipiv=i;
   for (int j=i;j<n;j++){
     if( abs(m[j][i])>piv){ 
       ipiv=j;
       piv= abs(m[j][i]);
       }
      }
//Permutar filas i y  ipiv 
     for (int j=i;j<n;j++) {
     double temp=m[i][j];
     m[i][j]=m[ipiv][j];
     m[ipiv][j]=temp;
     }  
     double temp = v[i];
     v[i]=v[ipiv];
     v[ipiv]=temp;

//Permutar vector permutacion  
     int itemp=per[i];           
     per[i]=per[ipiv];
     per[ipiv]=itemp;
       
//Eliminacion
  for (int j=i+1; j<n; j++){
    double   piv=m[j][i]/m[i][i]; //m[j][i] se anula despues de la primera elim.
    l[j][i]=piv;
    for (int k=i; k<n; k++){     
      m[j][k]=m[j][k]- piv*m[i][k];
      }
      v[j]=v[j]- piv*v[i];
    }
  }    
}

void backsus(int n,double**m, double* v, double* x){
//solucion de un sistema triangular superior por sustitucion hacia atras
for (int i=n-1;i>=0;i--){
   x[i]=v[i];
  for (int j= n-1;j>i;j--)
    x[i]=x[i]-m[i][j]*x[j];
    x[i]=x[i]/m[i][i];
  }
}     

void gauss_LU(int n, double** mat, double *b, double *x)
{
//n= dimension del sistema
//mat= matriz  coefficientes 
//b= vector terminos independientes

//Definir matriz L
double** l;
l = new double* [n];
for (int i=0;i<n;i++) l[i]=new double [n];


// Calcular sistema reducido y matriz l
gaussp(n,mat,b,l);

//Solucion del sistema por  substitucion hacia atras 
backsus(n,mat,b,x);

} 
*/
